In [1]:
%pylab
import numpy as np
import pylab as plb
import tifffile
In [2]:
#these functions are used for the image overlays.
def overlay_hue_saturation(overlay,bkgnd,
gama = 1.0,
color_map = plb.cm.jet,
alpha = 0.5):
from skimage import data, color, io, img_as_float
gama_correct = lambda x: x**gama
if len(np.shape(bkgnd)) == 2:
img_color = np.dstack((bkgnd, bkgnd, bkgnd))
else:
img_color = bkgnd
if len(np.shape(overlay)) == 2:
overlay_color = color_map(gama_correct(overlay))[:,:,:3]
mask = overlay>0
else:
overlay_color = overlay
mask = np.sum(overlay.astype(uint64),axis = 2)>0
img_hsv = color.rgb2hsv(img_color)
overlay_hsv = color.rgb2hsv(overlay_color)
img_hsv[..., 0] += mask*overlay_hsv[..., 0] * alpha
img_hsv[..., 1] += mask*overlay_hsv[..., 1] * alpha
img_overlay = color.hsv2rgb(img_hsv)
return(img_overlay)
def colorized_image(input_img,
map_point = 0.5,
color_map = plb.cm.jet,
alpha = 0.2,
gama = 1.0,
gain = 1.0):
from skimage import data, color, io, img_as_float
gama_correct = lambda x: x**gama
img = gain*input_img
img = gama_correct(img)
img[img>255] = 255
img = img.astype(uint8)
bg_color = np.dstack((img, img, img))
mask = img>0
map_img = mask*map_point
overlay_color = color_map(map_img)[:,:,:3]
img_hsv = color.rgb2hsv(bg_color)
overlay_hsv = color.rgb2hsv(overlay_color)
#make hue and saturation of img map to color
img_hsv[..., 0] += mask*overlay_hsv[..., 0] * alpha
img_hsv[..., 1] += mask*overlay_hsv[..., 1] * alpha
img_overlay = color.hsv2rgb(img_hsv)
return(img_overlay)
In [3]:
#load images
bkgnd_img = tifffile.TiffFile('./registration/65G06_brightfield_cuticle.tif').asarray()
b1_stack = tifffile.TiffFile('./registration/combined_stack_b1.tiff').asarray()
b2_stack = tifffile.TiffFile('./registration/combined_stack_b2.tiff').asarray()
b3_stack = tifffile.TiffFile('./registration/combined_stack_b3.tiff').asarray()
i1_stack = tifffile.TiffFile('./registration/combined_stack_i1.tiff').asarray()
i2_stack = tifffile.TiffFile('./registration/combined_stack_i2.tiff').asarray()
iii1_stack = tifffile.TiffFile('./registration/combined_stack_iii1.tiff').asarray()
iii3_stack = tifffile.TiffFile('./registration/combined_stack_iii3.tiff').asarray()
iii24_stack = tifffile.TiffFile('./registration/combined_stack_iii24.tiff').asarray()
hg1_stack = tifffile.TiffFile('./registration/combined_stack_hg1.tiff').asarray()
hg2_stack = tifffile.TiffFile('./registration/combined_stack_hg2.tiff').asarray()
hg3_stack = tifffile.TiffFile('./registration/combined_stack_hg3.tiff').asarray()
hg4_stack = tifffile.TiffFile('./registration/combined_stack_hg4.tiff').asarray()
tpd_stack = tifffile.TiffFile('./registration/combined_stack_tpd.tiff').asarray()
tpv_stack = tifffile.TiffFile('./registration/combined_stack_tpv.tiff').asarray()
pr_stack = tifffile.TiffFile('./registration/combined_stack_pr.tiff').asarray()
#make max projections
b1_img = np.max(b1_stack,axis = 0)
b2_img = np.max(b2_stack,axis = 0)
b3_img = np.max(b3_stack,axis = 0)
i1_img = np.max(i1_stack,axis = 0)
i2_img = np.max(i2_stack,axis = 0)
iii1_img = np.max(iii1_stack,axis = 0)
iii3_img = np.max(iii3_stack,axis = 0)
iii24_img = np.max(iii24_stack,axis = 0)
hg1_img = np.max(hg1_stack,axis = 0)
hg2_img = np.max(hg2_stack,axis = 0)
hg3_img = np.max(hg3_stack,axis = 0)
hg4_img = np.max(hg4_stack,axis = 0)
tpd_img = np.max(tpd_stack,axis = 0)
tpv_img = np.max(tpv_stack,axis = 0)
pr_img = np.max(pr_stack,axis = 0)
In [4]:
#colorize each muscle
b_map = plb.cm.Blues_r
i_map = plb.cm.Purples_r
iii_map = plb.cm.Reds_r
hg_map = plb.cm.Oranges_r
tp_map = plb.cm.Greens_r
pr_map = plb.cm.Greens_r
cvls = {'b1':0.01,'b2':0.2,'b3':0.4,
'i1':0.01,'i2':0.2,
'iii1':0.01,'iii3':0.2,'iii24':0.4,
'hg1':0.01,'hg2':0.2,'hg3':0.4,'hg4':0.6,
'tpd':0.01,'tpv':0.1,
'pr':0.01}
b1_colorized = colorized_image(b1_img,color_map = b_map,map_point = cvls['b1'],alpha = 1.0)
b2_colorized = colorized_image(b2_img,color_map = b_map,map_point = cvls['b2'],alpha = 1.0)
b3_colorized = colorized_image(b3_img,color_map = b_map,map_point = cvls['b3'],alpha = 1.0,gama = 0.65,gain = 40.0)
i1_colorized = colorized_image(i1_img,color_map = i_map,map_point = cvls['i1'],alpha = 1.0)
i2_colorized = colorized_image(i2_img,color_map = i_map,map_point = cvls['i2'],alpha = 1.0,gama = 0.7,gain = 10.0)
iii1_colorized = colorized_image(iii1_img,color_map = iii_map,map_point = cvls['iii1'],alpha = 1.0)
iii3_colorized = colorized_image(iii3_img,color_map = iii_map,map_point = cvls['iii3'],alpha = 1.0)
iii24_colorized = colorized_image(iii24_img,color_map = iii_map,map_point = cvls['iii24'],alpha = 1.0)
hg1_colorized = colorized_image(hg1_img,color_map = hg_map,map_point = cvls['hg1'],alpha = 1.0,gain = 2.0)
hg2_colorized = colorized_image(hg2_img,color_map = hg_map,map_point = cvls['hg2'],alpha = 1.0)
hg3_colorized = colorized_image(hg3_img,color_map = hg_map,map_point = cvls['hg3'],alpha = 1.0)
hg4_colorized = colorized_image(hg4_img,color_map = hg_map,map_point = cvls['hg4'],alpha = 1.0)
tpd_colorized = colorized_image(tpd_img,color_map = tp_map,map_point = cvls['tpd'],alpha = 1.0)
tpv_colorized = colorized_image(tpv_img,color_map = tp_map,map_point = cvls['tpv'],alpha = 1.0)
pr_colorized = colorized_image(pr_img,color_map = pr_map,map_point = cvls['pr'],alpha = 1.0)
#combine the images using addition
img_sum = (tpd_colorized+
tpv_colorized+
b1_colorized+
b2_colorized+
b3_colorized+
i1_colorized+
i2_colorized+
iii1_colorized+
iii3_colorized+
iii24_colorized+
hg1_colorized+
hg2_colorized+
hg3_colorized+
hg4_colorized+
pr_colorized)*255
img_sum[img_sum>255] = 255
img_sum = img_sum.astype(uint8)
imshow(overlay_hue_saturation(img_sum[:,:,:3],bkgnd_img,alpha = 1.0))
display(gcf());close()
I recursively apply the function 'mask_stack' to a list of colorized muscle stacks to create a visualization of the muscle population. The order of muscles in the list is important since this determines the way the muscles overlap in the final image. The colorized image of the muscles are then mixed into the greyscale image of the cuticle.
In [5]:
def mask_stack(top_img,bottom_img,alpha = 0.5,gain = 15):
from skimage import data, color, io, img_as_float
top_hsv = color.rgb2hsv(top_img)
mask = top_hsv[...,2]/np.max(top_hsv[...,2])
mask *= gain
mask[mask>1] = 1
mask[mask<1] = 0
mask = np.dstack((mask,mask,mask))
inverse = 1-mask*alpha
return top_img*mask + bottom_img*inverse
stacked_muscles = reduce(mask_stack,
[tpd_colorized,
tpv_colorized,
pr_colorized,
b2_colorized,
i1_colorized,
hg4_colorized,
iii24_colorized,
i2_colorized,
hg2_colorized,
hg3_colorized,
iii3_colorized,
iii1_colorized,
b1_colorized,
b3_colorized,
hg1_colorized])
gama = 1.5
gama_correct = lambda x: x**gama
cut_gamma = gama_correct(bkgnd_img)+600
cuticle = (np.dstack((cut_gamma,cut_gamma,cut_gamma))*0.03).astype(uint8)
muscles = (stacked_muscles*0.9*255)
sumimg = cuticle+muscles
sumimg[sumimg>255] = 255
sumimg = sumimg.astype(uint8)
imshow(sumimg)
display(gcf());close()
In [6]:
tifffile.imsave('stacked_muscles.tiff',uint8(stacked_muscles*255))
tifffile.imsave('muscle_viz.tiff',sumimg)
To define the regions of interest I find the 2d contours around each muscle. For simplicity I then resample the contours to 20 points: hopefully this should allow us to average in the map extracted from other confocal stacks in the future. We then save the lists of contours to use in the construction of the model.
In [7]:
mimgs = {'b1':b1_img,
'b2':b2_img,
'b3':b3_img,
'i1':i1_img,
'i2':i2_img,
'iii1':iii1_img,
'iii3':iii3_img,
'iii4':iii24_img,
'hg1':hg1_img,
'hg2':hg2_img,
'hg3':hg3_img,
'hg4':hg4_img,
'tpv':tpv_img,
'tpd':tpd_img,
'pr':pr_img}
model_data = dict()
#construct model
from skimage import measure
from scipy.interpolate import griddata
figure(figsize(5,5))
imshow(sumimg)
for mkey in mimgs.keys():
img = mimgs[mkey]
contours = [measure.find_contours(img,1)[0]]
for n, contour in enumerate(contours):
clen = len(contour[:,0])
x = griddata(np.arange(clen),contour[:,0],np.linspace(clen,20))
x = x[~isnan(x)]
x = hstack((x,x[0]))
y = griddata(np.arange(clen),contour[:,1],np.linspace(clen,20))
y = y[~isnan(y)]
y = hstack((y,y[0]))
model_data[mkey] = vstack((y,x))
plot(y, x,linewidth=0.5,color = 'w')
gca().set_xbound((0,1024))
gca().set_ybound((0,1024))
display(gcf());close()
In [8]:
class Basis(dict):
def __setitem__(self,key,item):
"""overload the __setitem__ method of dict so the transform and inverse
transform matrices are computed when the basis vectors are changed"""
try:
if key in ['a1','a2']:
dict.__setitem__(self,key,item)
A = np.vstack((self['a1'],self['a2'])).T
A_inv = np.linalg.inv(A)
self['A'] = A
self['A_inv'] = A_inv
else:
dict.__setitem__(self,key,item)
except KeyError:
dict.__setitem__(self,key,item)
class GeometricModel(object):
def __init__(self,lines,basis):
self.lines = lines
self.basis = basis
## put lines in barycentric coords
self.barycentric = dict()
for key in self.lines.keys():
coords = dot(self.basis['A_inv'],(self.lines[key]-self.basis['p'][:,newaxis]))
self.barycentric[key] = coords.T
def coords_from_basis(self,basis):
ret = dict()
for key in self.barycentric.keys():
coords = np.dot(basis['A'],(self.barycentric[key]).T)+basis['p'][:,newaxis]
ret[key] = coords
return(ret)
class ModelView(object):
def __init__(self,model):
self.model = model
def plot(self,basis,**kwargs):
lines = self.model.coords_from_basis(basis)
plot_args = {}
plot_args['plot_frame'] = kwargs.pop('plot_frame',True)
plot_args['frame_head_width'] = kwargs.pop('frame_head_width',20)
plot_args['contour_color'] = kwargs.pop('contour_color','w')
kwargs['color'] = plot_args['contour_color']
for line in lines.values():
plot(line[0,:],line[1,:], **kwargs)
if plot_args['plot_frame']:
p = basis['p']
a1 = basis['a1']
a2 = basis['a2']
kwargs['color'] = 'g'
kwargs['head_width'] = plot_args['frame_head_width']
arrow(p[0],p[1],a1[0],a1[1],**kwargs)
kwargs['color'] = 'b'
kwargs['head_width'] = plot_args['frame_head_width']
arrow(p[0],p[1],a2[0],a2[1],**kwargs)
Now I construct a model and plot it in the reference frame of the confocal image: the position of the three bristles that form the reference frame are hard coded here.
In [9]:
figure(figsize = (10,10))
imfile = tifffile.TiffFile('im_overlay.tiff')
sumimg = imfile.asarray()
###add position of large setae
model_data['e1'] = np.array([[ 170.02788104, 326.71685254],
[ 380.98203222, 919.92627014]])
model_data['e2'] = array([[ 172.83333333, 332.83333333],
[ 551.5 , 164.83333333]])
e1 = model_data['e1']
e2 = model_data['e2']
muscles = dict()
for key in model_data.keys():
if not(key in ['e1','e2']):
muscles[key] = model_data[key]
confocal_basis = Basis()
confocal_basis['a1'] = e2[1]-e2[0]
confocal_basis['a2'] = e1[1]-e2[0]
confocal_basis['p'] = e2[0]
thorax = GeometricModel(muscles,confocal_basis)
thorax_view = ModelView(thorax)
imshow(sumimg)
import copy
thorax_view.plot(thorax.basis)
gca().axis('tight')
display(gcf());close()
Using the imageAnalysis.py program I manually fit the reference vectors to the gcamp imaging stacks. With the extracted point corespondences, the images from the gcamp movies can then be warped back into the same frame as the confocal stack. Warp and save will perform this on a given fly. I then plot the max projection images from a squadron of eight flies.
In [10]:
# define boolean masks for each muscle
b1_mask = b1_img >0
b2_mask = b2_img >0
b3_mask = b3_img >0
i1_mask = i1_img >0
i2_mask = i2_img >0
iii1_mask = iii1_img >0
iii3_mask = iii3_img >0
iii24_mask = iii24_img >0
hg1_mask = hg1_img >0
hg2_mask = hg2_img >0
hg3_mask = hg3_img >0
hg4_mask = hg4_img >0
tpv_mask = tpv_img >0
tpd_mask = tpd_img >0
pr_mask = pr_img > 0
#Masks for each line
mU56_X_J133 = b1_mask | b2_mask | b3_mask \
| i1_mask | i2_mask \
| iii1_mask | iii24_mask |iii3_mask \
| hg1_mask | hg2_mask | hg3_mask | hg4_mask \
| tpd_mask | tpv_mask \
| pr_mask
mU82_X_J18 = b1_mask | b2_mask | b3_mask \
| i1_mask | i2_mask \
| iii1_mask |iii3_mask \
| hg2_mask | hg3_mask | hg4_mask \
mU82_X_J100 = iii24_mask | hg2_mask
mU82_X_J22 = b1_mask | b2_mask | b3_mask \
|iii3_mask \
| hg1_mask | hg2_mask | hg3_mask | hg4_mask \
mU82_X_J23 = b1_mask | b2_mask | b3_mask \
| iii1_mask |iii3_mask \
| hg2_mask | hg4_mask \
mU82_X_J121 = tpd_mask | tpv_mask
In [11]:
def blur_by_depth(stack):
from scipy.ndimage.filters import gaussian_filter
wf = lambda x: (1.0/(x+30))+1.0
return np.array([gaussian_filter(stack[i],i/2.0 + 3.0)*wf(i) for i in range(shape(stack)[0])])
def model_muscle(stack):
return sum(blur_by_depth((stack>0).astype(float)),axis = 0)
b1_model = model_muscle(b1_stack)
b2_model = model_muscle(b2_stack)
b3_model = model_muscle(b3_stack)
i1_model = model_muscle(i1_stack)
i2_model = model_muscle(i2_stack)
iii1_model = model_muscle(iii1_stack)
iii3_model = model_muscle(iii3_stack)
iii24_model = model_muscle(iii24_stack)
hg1_model = model_muscle(hg1_stack)
hg2_model = model_muscle(hg2_stack)
hg3_model = model_muscle(hg3_stack)
hg4_model = model_muscle(hg4_stack)
tpd_model = model_muscle(tpd_stack)
tpv_model = model_muscle(tpv_stack)
pr_model = model_muscle(pr_stack)
\(X\) is an matrix of \(m\) pixels by \(n\) frames constructed from the signals from \(k\) muscles so that: \(X = WB\) where \(W\) is a \(m\) by \(k\) mixing matrix and \(B\) is a \(k\) by \(n\) matrix of muscle fluorescence over time. We find the best fit to \(B\) as \(W^\dagger X = B\).
In [12]:
import cv2
from scipy.ndimage.morphology import binary_dilation
bk = binary_dilation(mU56_X_J133,iterations = 60).astype(float)
from scipy.ndimage.filters import gaussian_filter
bk = gaussian_filter(bk,50)
mlist = [b1_model,b2_model,b3_model,
i1_model,i2_model,
iii1_model,iii3_model,iii24_model,
hg1_model,hg2_model,hg3_model,hg4_model,
tpd_model,tpv_model,pr_model]#,bk]
model_mtrx = [cv2.pyrDown(cv2.pyrDown(ar)) for ar in mlist]
W = np.array([m.ravel() for m in model_mtrx]).T
In [13]:
W_ = W#= W/np.max(W,axis = 0)
imshow(sum(W_.reshape(256,256,-1),axis = 2),cmap = cm.gray)
display(gcf());close()
In [14]:
#Groups for screen
import flylib as flb
import db_access as dba
import cv2
import cPickle
fly_db = dba.get_db()
U82_X_J121 = flb.Squadron(fly_db,[271,272,273,274,275,276])
U82_X_J18 = flb.Squadron(fly_db,[265,267,268,269,270])##266
U82_X_J100 = flb.Squadron(fly_db,[259,260,261,262,263,264])
U82_X_J22 = flb.Squadron(fly_db,[255,256,257,258])
U82_X_J23 = flb.Squadron(fly_db,[278,279,280,281,282]) ##277
U56_X_J133 = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])
swarm_dict = {'U56_X_J133':U56_X_J133,
'U82_X_J18':U82_X_J18,
'U82_X_J100':U82_X_J100,
'U82_X_J22':U82_X_J22,
'U82_X_J23':U82_X_J23,
'U82_X_J121':U82_X_J121}
import h5py
testfly = swarm_dict['U56_X_J133'].flies[4]
f = h5py.File(testfly.fly_path + "warped_imgs.hdf5", "r")
X_warped = np.array(f['warped_imgs'])
X = X_warped.T.astype(float)
print shape(X)
print shape(W_)
print numpy.linalg.matrix_rank(W_)
W_pin = numpy.linalg.pinv(W_)
B = dot(W_pin,X)
In [17]:
lns = plot(B.T)
gca().set_xbound(500,2500)
display(gcf());close()
In [20]:
mlist = [b1_model,b2_model,b3_model,
i1_model,i2_model,
iii1_model,iii3_model,iii24_model,
hg1_model,hg2_model,hg3_model,hg4_model,
tpd_model,tpv_model,pr_model,bk]
model_mtrx = [cv2.pyrDown(cv2.pyrDown(ar)) for ar in mlist]
W = np.array([m.ravel() for m in model_mtrx]).T
W_ = W
W_pin = numpy.linalg.pinv(W_)
B = dot(W_pin,X)
In [25]:
lns = plot(B.T)
gca().set_xbound(500,2500)
gca().set_ybound(-0.5,1.5)
display(gcf());close()
In []:
In []:
In []:
In []:
In []:
In []:
In []:
In [16]:
X_est = dot(W_,B[:,:5000])
#res = X[:,:5000] - X_est
savetiff = hstack((transpose(X_warped[:5000,:].reshape(5000,256,256),(0,2,1))*5,
transpose(X_est.T.reshape(5000,256,256),(0,2,1))*5))
tifffile.imsave('model_fit4.tiff',transpose(savetiff.astype(uint8),(0,2,1)))
In []:
In []:
In []:
In []:
In []:
In []:
In []:
In [16]:
savetiff = hstack((transpose(X_warped[:5000,:].reshape(5000,256,256),(0,2,1))*5,
transpose(X_est.T.reshape(5000,256,256),(0,2,1))*5))
#savetiff += np.min(savetiff)
#savetiff[savetiff > 256] = 256
tifffile.imsave('model_fit5.tiff',transpose(savetiff.astype(uint8),(0,2,1)))
In []:
In []:
In [61]:
hist(X_est.T[:5000,:].ravel())
Out[61]:
In []:
max(savetiff.ravel())
In [86]:
for i in range(16):
subplot(16,1,i+1,sharex = gca(),sharey = gca())
plot(B[i,:])
In [51]:
res = X[:,:5000] - X_est
In [56]:
imshow(res[:,0].reshape(256,256),cmap = cm.gray)
Out[56]:
In []:
In []:
In []:
In []:
In []:
In [10]:
import db_access as dba
import flylib as flb
import cv2
import cPickle
output_shape = (1024,1024)
figure(figsize = (8,16))
fly_db = dba.get_db()
swarm = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])
for i,fly in enumerate(swarm.flies):
#fly = swarm.flies[0]
tiffdata = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
pkname = fly.fly_path + '/basis_fits.cpkl'
#print 'here'
f = open(pkname)
img_data = cPickle.load(f)
f.close()
test_basis = Basis()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
maximg = np.max(tiffdata,axis=0)
col = mod(i,2)
row = i/2
subplot2grid((4,2),(row,col))
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
imshow(cv2.warpAffine(maximg.T,A,output_shape),cmap = cm.gray)
thorax_view.plot(thorax.basis,alpha = 0.5)
gca().set_ybound(0,1024)
gca().set_xbound(0,1024)
gca().set_xticks([])
gca().set_yticks([])
display(gcf());close()
Now I start with a preliminary stab at group analysis.
In [11]:
#Groups for screen
U82_X_J121 = flb.Squadron(fly_db,[271,272,273,274,275,276])
U82_X_J18 = flb.Squadron(fly_db,[265,267,268,269,270])##266
U82_X_J100 = flb.Squadron(fly_db,[259,260,261,262,263,264])
U82_X_J22 = flb.Squadron(fly_db,[255,256,257,258])
U82_X_J23 = flb.Squadron(fly_db,[278,279,280,281,282]) ##277
U56_X_J133 = flb.Squadron(fly_db,[244,243,242,241,240,239,238,237])
# define boolean masks for each muscle
b1_mask = b1_img >0
b2_mask = b2_img >0
b3_mask = b3_img >0
i1_mask = i1_img >0
i2_mask = i2_img >0
iii1_mask = iii1_img >0
iii3_mask = iii3_img >0
iii24_mask = iii24_img >0
hg1_mask = hg1_img >0
hg2_mask = hg2_img >0
hg3_mask = hg3_img >0
hg4_mask = hg4_img >0
tpv_mask = tpv_img >0
tpd_mask = tpd_img >0
#Masks for each line
mU56_X_J133 = b1_mask | b2_mask | b3_mask \
| i1_mask | i2_mask \
| iii1_mask | iii24_mask |iii3_mask \
| hg1_mask | hg2_mask | hg3_mask | hg4_mask \
| tpd_mask | tpv_mask
mU82_X_J18 = b1_mask | b2_mask | b3_mask \
| i1_mask | i2_mask \
| iii1_mask |iii3_mask \
| hg2_mask | hg3_mask | hg4_mask \
mU82_X_J100 = iii24_mask | hg2_mask
mU82_X_J22 = b1_mask | b2_mask | b3_mask \
|iii3_mask \
| hg1_mask | hg2_mask | hg3_mask | hg4_mask \
mU82_X_J23 = b1_mask | b2_mask | b3_mask \
| iii1_mask |iii3_mask \
| hg2_mask | hg4_mask \
mU82_X_J121 = tpd_mask | tpv_mask
#for plotting
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
In [12]:
i1_select = i1_mask & ~(b2_mask | b3_mask | b1_mask)
i2_select = i2_mask & ~(iii3_mask | hg4_mask | hg1_mask)
b1_select = b1_mask & ~ (i1_mask | b2_mask | b3_mask)
b2_select = b2_mask & ~ (i1_mask | b1_mask | b3_mask | tpd_mask)
b3_select = b3_mask & ~ (i1_mask | b2_mask | b1_mask | tpd_mask)
iii1_select = iii1_mask & ~ (b2_mask | tpd_mask | tpv_mask | iii24_mask)
iii3_select = iii3_mask & ~ (hg1_mask | hg4_mask | hg3_mask | i2_mask)
iii23_select = iii24_mask & ~ (hg1_mask | hg4_mask | hg3_mask | i2_mask | iii3_mask)
hg1_select = hg1_mask & ~ (hg4_mask | hg3_mask | iii3_mask | iii24_mask)
hg2_select = hg2_mask & ~ (hg4_mask)
hg3_select = hg3_mask & ~ (hg2_mask)
hg4_select = hg4_mask & ~ (iii3_mask | iii24_mask)
tp_select = (tpd_mask | tpv_mask) & ~ (i1_mask | b3_mask | iii1_mask | i2_mask)
figure(figsize(12,4))
for i,mask in enumerate([b1_select,b2_select,b3_select,
i1_select,i2_select,
iii1_select, iii3_select,iii23_select,
hg1_select,hg2_select,hg3_select,hg4_select,
tp_select]):
subplot(2,7,i+1)
imshow(mask)
thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False,contour_color = 'white',alpha = 0.5)
gca().set_xticklabels([]);gca().set_yticklabels([]);
gca().set_xbound((0,1024));gca().set_ybound((0,1024))
display(gcf());close()
In [19]:
figure()
imshow((hg4_model)+b1_model,cmap=cm.gray)
Out[19]:
In [21]:
imshow(sum(model_mtrx,axis = 0),cmap = cm.gray)
Out[21]:
In [24]:
plot(B.T)
Out[24]:
In [98]:
imshow(X_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(sum(res.T.reshape(-1,256,256),axis = 0),cmap = cm.gray)
Out[98]:
In []:
In [14]:
figure(figsize(8,10))
fly = U56_X_J133.flies[1]
expmnt = fly.experiments.values()[0]
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])
rwing = rad2deg(np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])/5.0)
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
imgs = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()
ax = subplot2grid((14,6),(0,0),colspan = 5)
plot(times[:30000],ypos[:30000],rasterized=True)
locator_params(axis = 'y', nbins = 3)
subplot2grid((14,6),(1,0),colspan = 5,sharex = ax)
plot(times[:30000],rwing[:30000],rasterized=True)
locator_params(axis = 'y', nbins = 3)
gca().set_ybound(45,90)
for i,mask in enumerate([b1_select,b2_select,b3_select,
i1_select,i2_select,
iii1_select,iii3_select,iii23_select,
hg1_select,hg2_select,hg3_select,hg4_select]):
subplot2grid((14,6),(2+i,0),colspan = 5,sharex = ax)
mmask = cv2.pyrDown(cv2.pyrDown((mask).astype(float)))>0
plot(times[:30000],mean(imgs[:,mmask],axis = 1),rasterized=True)
if i<11:
gca().set_xticklabels([])
locator_params(axis = 'y', nbins = 3)
subplot2grid((14,6),(2+i,5))
imshow(mask)
thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False,contour_color = 'white',alpha = 0.5)
gca().set_xticklabels([]);gca().set_yticklabels([]);
gca().set_xbound((0,1024));gca().set_ybound((0,1024))
subplots_adjust(hspace = 0.1)
subplots_adjust(bottom = 0.05)
subplots_adjust(top = 0.99)
subplots_adjust(right = 0.99)
subplots_adjust(left = 0.05)
subplots_adjust(wspace = 0.00)
ax.set_xbound(10,140)
display(gcf());close()
In [90]:
In []:
In []:
In [97]:
from scipy.optimize import nnls
blist = list()
for i in range(1000):
blist.append(nnls(W,X[:,i]))
In [98]:
B_nn = array([x[0] for x in blist])
In [99]:
X_nn_est = dot(W,B_nn.T)
res_nn = X[:,:1000] - X_nn_est
res = X[:,:1000] - X_est
In [100]:
imshow(X_nn_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(sum(res_nn.T.reshape(-1,256,256),axis = 0),cmap = cm.gray)
figure()
imshow(X_est.T.reshape(-1,256,256)[0],cmap = cm.gray)
figure()
imshow(res.T.reshape(-1,256,256)[0],cmap = cm.gray)
Out[100]:
In [104]:
plot(B.T)
Out[104]:
In [103]:
close()
In []:
In []:
In []:
In []:
In [94]:
#this block will load the metadata and warp the images/movies into a common reference frame.
swarm_dict = {'U56_X_J133':U56_X_J133,
'U82_X_J18':U82_X_J18,
'U82_X_J100':U82_X_J100,
'U82_X_J22':U82_X_J22,
'U82_X_J23':U82_X_J23,
'U82_X_J121':U82_X_J121}
def warp_fly_movie(fly):
fname = fly.fly_path + 'meta_data.cpkl'
f = open(fname,'rb');data = cPickle.load(f);f.close()
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname);img_data = cPickle.load(f);f.close()
test_basis = Basis()
#print data.keys()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
#compose back to Ca image size
A = vstack((A,[0,0,1]))
S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
Ap = dot(S,A)[:2,:3]
testmovie = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
X_warped = [cv2.warpAffine(testmovie[frame].T,Ap,(256,256)).ravel() for frame in range(shape(testmovie)[0])]
X_warped = np.array(X_warped)
import h5py
f = h5py.File(fly.fly_path + "warped_imgs.hdf5", "a")
f['warped_imgs'] = X_warped
f.flush()
f.close()
In [95]:
for swarm in swarm_dict.values():
for fly in swarm.flies:
# os.remove(fly.fly_path + "warped_imgs.hdf5")
print fly.fly_num
warp_fly_movie(fly)
In [87]:
plot(B[:,:5000].T)
Out[87]:
In [88]:
X_est = dot(W,B[:,:5000])
In [58]:
print mean(B[:,:5000],axis = 1)
In [60]:
from scipy.optimize import nnls
In [133]:
Out[133]:
In [56]:
imshow(Out[51][0])
Out[56]:
In [15]:
#this block will load the metadata and warp the images/movies into a common reference frame.
swarm_dict = {'U56_X_J133':U56_X_J133,
'U82_X_J18':U82_X_J18,
'U82_X_J100':U82_X_J100,
'U82_X_J22':U82_X_J22,
'U82_X_J23':U82_X_J23,
'U82_X_J121':U82_X_J121}
img_dict = dict()
output_shape = (256,256)
import cPickle
for swarm_key in swarm_dict.keys():
swarm = swarm_dict[swarm_key]
corr_list = list()
dot_list = list()
max_list = list()
wing_list = list()
norm_list = list()
mean_list = list()
movie_list = list()
rwing_list = list()
ypos_list = list()
for fly in swarm.flies:
try:
#print(fly.fly_num)
fname = fly.fly_path + 'meta_data.cpkl'
#f = open(fname,'rb');data = cPickle.load(f)['cl_stim_metadata'];f.close()
f = open(fname,'rb');data = cPickle.load(f);f.close()
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname);img_data = cPickle.load(f);f.close()
test_basis = Basis()
#print data.keys()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
#compose back to Ca image size
A = vstack((A,[0,0,1]))
S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
Ap = dot(S,A)[:2,:3]
#warp the images
dot_list.append(cv2.warpAffine(data['sin_stim_metadata']['dot_img'].T,Ap,output_shape))
corr_list.append(cv2.warpAffine(data['sin_stim_metadata']['corr_img'].T,Ap,output_shape))
max_list.append(cv2.warpAffine(data['sin_stim_metadata']['max_img'].T,Ap,output_shape))
mean_list.append(cv2.warpAffine(data['sin_stim_metadata']['mean_img'].T,Ap,output_shape))
norm_list.append(cv2.warpAffine(data['sin_stim_metadata']['norm_img'].T,Ap,output_shape))
#now warp the movies
new_shape = (shape(data['ave_movie'])[0],output_shape[0],output_shape[1])
warped_ave_movie = zeros(new_shape)
for frame in arange(shape(warped_ave_movie)[0]):
warped_ave_movie[frame,:,:] = cv2.warpAffine(data['ave_movie'][frame,:,:].T,Ap,output_shape)
movie_list.append(warped_ave_movie/mean_list[-1])
rwing_list.append(data['ave_rwing'])
ypos_list.append(data['ave_ypos'])
except Exception as inst:
print inst
img_dict[swarm_key] = {'dot_list':dot_list,
'corr_list':corr_list,
'mean_list':mean_list,
'norm_list':norm_list,
'max_list':max_list,
'movie_list':movie_list,
'rwing_list':rwing_list,
'ypos_list':ypos_list}
In [16]:
ndat = np.array(img_dict[swarm_key]['norm_list'][:]).ravel()
ddat = np.array(img_dict[swarm_key]['dot_list'][:]).ravel()
cdat = np.array(img_dict[swarm_key]['corr_list'][:]).ravel()
norm_max = np.percentile(ndat,95)
dot_max = np.percentile(ddat,95)
corr_max = np.percentile(cdat,95)
norm_min = np.percentile(ndat,0)
dot_min = np.percentile(ddat,5)
corr_min = np.percentile(cdat,5)
def summarize_swarm(swarm_key,swarm_mask,axlist):
axes(axlist[0])
imshow(swarm_mask)
thorax_view.plot(thorax.basis,lw = 1.0,plot_frame = False)
axes(axlist[1])
imshow(np.mean(img_dict[swarm_key]['norm_list'],axis = 0),cmap = cm.gray,vmin = 0, vmax = norm_max)
thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)
axes(axlist[2])
imshow(np.mean(img_dict[swarm_key]['dot_list'],axis = 0),cmap = cm.gray,vmin = dot_min, vmax = dot_max)
thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)
axes(axlist[3])
imshow(np.mean(img_dict[swarm_key]['corr_list'],axis = 0),cmap = cm.gray,vmin = corr_min, vmax = corr_max)
thorax_view.plot(little_frame,lw = 1.0,plot_frame = False,contour_color = 'orange',alpha = 0.5)
figure(figsize(8,10))
axlist = [subplot(5,4,i) for i in range(1,5)]
summarize_swarm('U56_X_J133',mU56_X_J133,axlist)
axlist[0].set_ylabel("U56 X J133")
#axlist[1].set_title(r"$\sqrt{diag(XX^\top)}$")
#axlist[2].set_title(r"$Xs^\top$")
#axlist[3].set_title(r"$\frac{Xs^\top}{\sqrt{diag(XX^\top)}||s||}$",y= 1.08)
axlist[1].set_title(r"pixel std")
axlist[2].set_title(r"dot image")
axlist[3].set_title(r"correlation image")
axlist = [subplot(5,4,i) for i in range(5,9)]
summarize_swarm('U82_X_J18',mU82_X_J18,axlist)
axlist[0].set_ylabel("U82 X J18")
axlist = [subplot(5,4,i) for i in range(9,13)]
summarize_swarm('U82_X_J100', mU82_X_J100,axlist)
axlist[0].set_ylabel("U82 X J100")
axlist = [subplot(5,4,i) for i in range(13,17)]
summarize_swarm('U82_X_J22', mU82_X_J22,axlist)
axlist[0].set_ylabel("U82 X J22")
axlist = [subplot(5,4,i) for i in range(17,21)]
summarize_swarm('U82_X_J23', mU82_X_J23,axlist)
axlist[0].set_ylabel("U82 X J23")
for ax in gcf().axes:
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_xticks([]); ax.set_yticks([])
for i in [1,2,3,
5,6,7,
9,10,11,
13,14,15,
17,18,19]:
ax = gcf().axes[i]
ax.set_xbound(15,190);ax.set_ybound(60,236)
for i in [0,4,8,12,16]:
ax = gcf().axes[i]
ax.set_xbound(15*4,190*4);ax.set_ybound(60*4,235*4)
#tight_layout()
display(gcf());#close()
In [17]:
swarm_dict = {'U56_X_J133':U56_X_J133,
'U82_X_J18':U82_X_J18,
'U82_X_J100':U82_X_J100,
'U82_X_J22':U82_X_J22,
'U82_X_J23':U82_X_J23,
'U82_X_J121':U82_X_J121}
m_dict = {'U56_X_J133':mU56_X_J133,
'U82_X_J18':mU82_X_J18,
'U82_X_J100':mU82_X_J100,
'U82_X_J22':mU82_X_J22,
'U82_X_J23':mU82_X_J23,
'U82_X_J121':mU82_X_J121}
def plot_group_sig(sig):
sig_mean = np.mean(sig,axis = 0)
sig_std = np.std(sig,axis = 0)
times = np.linspace(0,2,60*2)
#fill_between(times,sig_mean+sig_std,sig_mean-sig_std,alpha = 0.5)
plot(times,sig.T,color = 'k',alpha = 0.5)#sig_mean)
plot(times,sig_mean,lw = 2.0,alpha = 0.8)
figure(figsize = (8,10))
for j,swarm_key in enumerate(swarm_dict.keys()):
#swarm_key = 'U82_X_J23'
subplot2grid((13,7),(0,j),colspan = 1)
imshow(m_dict[swarm_key])
thorax_view.plot(thorax.basis,lw = 0.5,alpha = 0.5,plot_frame = False)
gca().set_xbound(0,1024);gca().set_ybound(0,1024)
wing_angle = rad2deg(np.array(img_dict[swarm_key]['rwing_list'])/5)
subplot2grid((13,7),(1,j),colspan = 1)
plot_group_sig(wing_angle)
gca().set_ybound(45,75)
def plot_roi(mask):
mmask = cv2.pyrDown(cv2.pyrDown(mask.astype(float)))>0
mdata = np.mean(np.array(img_dict[swarm_key]['movie_list'])[:,:,mmask],axis=2)
plot_group_sig(mdata)
gca().set_ybound(0.9,1.1)
def imshow_roi(mask):
mmask = cv2.pyrDown(cv2.pyrDown(mask.astype(float)))>0
imshow(mmask)
thorax_view.plot(little_frame,lw = 0.5,plot_frame = False,contour_color = 'white',alpha = 0.5)
gca().set_xticklabels([]);gca().set_yticklabels([]);
gca().set_xbound(15,190);gca().set_ybound(60,236)
for i,mask in enumerate([b1_select,
b2_select,
b3_select,
i1_select,
i2_select,
iii1_select,
iii3_select,
hg1_select,
hg2_select,
hg3_select,
hg4_select,
tp_select]):
subplot2grid((14,7),(i+2,j),colspan = 1)
plot_roi(mask)
subplot2grid((14,7),(i+2,6))
imshow_roi(mask)
for ax in gcf().axes:
ax.set_yticklabels([])
ax.set_xticklabels([])
subplots_adjust(hspace = 0.1)
subplots_adjust(bottom = 0.05)
subplots_adjust(top = 0.99)
subplots_adjust(right = 0.99)
subplots_adjust(left = 0.05)
subplots_adjust(wspace = 0.00)
display(gcf())
In [20]:
print(shape(mask_mtrx))
print(numpy.linalg.matrix_rank(mask_mtrx))
In [21]:
imshow(sum(mask_mtrx,axis = 0).reshape(256,256))
Out[21]:
In [24]:
In [42]:
X_mean = np.mean(X_warped.astype(float),axis = 0)
#X_Xbar = X_warped/X_mean
X_norm = np.linalg.norm(X_warped.astype(float),axis = 0)
In [59]:
X_Xnorm = (X_warped.astype(float)-X_mean)[~(xwarped == 0)]/X_norm
In [58]:
hist(X_Xnorm[~isnan(X_Xnorm)].ravel())
Out[58]:
In [60]:
del(np)
In [47]:
imshow(X_mean.reshape(256,256))
Out[47]:
In [213]:
gc.collect()
Out[213]:
In [200]:
plot(X_Xbar[500])
Out[200]:
In [129]:
imshow(fit_movie.T[9].reshape(256,256),cmap = cm.gray)
Out[129]:
In [134]:
tifffile.imsave('fittest2.tiff',fit_movie.T.reshape(10000,256,256).astype(uint8)*10)
In [138]:
tifffile.imsave('fittest_X.tiff',X_warped.reshape(10000,256,256).astype(uint8))
In []:
In [85]:
#functions to summarize the sin stimulus protocol
def collect_cycles(ypos,stim_cond,trial_selection_mask,flight_mask):
cycles = list()
trial_selector = zeros_like(ypos).astype(bool)
cycle_selection_mask = ypos>1.5
#trial_selection_mask = around(stim_cond)==2
trial_idxs = flb.idx_by_thresh(trial_selection_mask,0.5)
for trial_idx in trial_idxs:
if not(np.sum((~flight_mask)[trial_idx])):
trial_selector[:] = 0
trial_selector[trial_idx] = 1
cycle_idxs = flb.idx_by_thresh(trial_selector & cycle_selection_mask)
cycles.extend([arange(cycle_idxs[i][0],cycle_idxs[i+1][0])
for i in range(len(cycle_idxs)-1)])
return cycles
def cycle_average(sig,cycle_idxs,times,flight_mask):
"""function extracts the cycle average of the signal given a list of idxs where the cycles
happen, and a mask to exclude non-flight bouts. Resamples so that cycle samples line up.
returns a dictionary of summary stats"""
trial_idxs = [np.arange(x[0],x[0]+210) for x in cycle_idxs]# if (len(x)>0 and (sum(~flight_mask[x[0]:x[0]+200])==0))]
trial_times = [times[idx]-times[idx[0]] for idx in trial_idxs if sum(~flight_mask[idx])==0]
trial_data = [sig[idx] for idx in trial_idxs if sum(~flight_mask[idx])==0]
import scipy.interpolate
fun = scipy.interpolate.griddata
time_resamp = linspace(0,2.,60*2.)
data_resamp = list()
for td,tt in zip(trial_data,trial_times):
f = scipy.interpolate.interp1d(tt, td, kind='linear', axis=0)
data_resamp.append(f(time_resamp))
retdict = {'time_resamp':time_resamp,
'sig_resamp':data_resamp}
return retdict
def wing_correlations(sig,imgs,fly = None):
#calculate some summary images including the correlation image
mean_img = mean(imgs,axis=0)
zeroed_img = imgs - mean_img
max_img = np.max(imgs,axis = 0)
norm_img = numpy.linalg.norm(zeroed_img,axis = 0)
norm_sig = numpy.linalg.norm(sig-np.mean(sig))
dot_img = sum(zeroed_img*(sig[:,newaxis,newaxis]-np.mean(sig)),axis = 0)
corr_img = dot_img/(norm_img*norm_sig)
return {'mean_img':mean_img,
'max_img':max_img,
'norm_img':norm_img,
'dot_img':dot_img,
'corr_img':corr_img}
def calc_fly_stats(fly):
#only did one experiment on these flies
expmnt = fly.experiments.values()[0]
nframes = 30000
#load the movie
movie = np.array(expmnt.exp_record['tiff_data']['images'][:nframes])
flight_mask = np.array(expmnt.exp_record['flight_mask'][:nframes])
#Map some variables
stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])[:nframes]
lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])[:nframes]
rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])[:nframes]
lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])[:nframes]
aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])[:nframes]
xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])[:nframes]
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])[:nframes]
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])[:nframes]
#idxs associated with sin type stimulus
trial_selection_mask = stim_cond>0
cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask,flight_mask)
#concatinate the stimuli idxs
sin_idxs = hstack(cycle_idxs)
#select the frames from the imaging stack
#that correspond to the sin stimulus
imgs = movie[sin_idxs,:,:]
sig = rwing[sin_idxs]
#correlate with the right wingstroke data
sin_stim_metadata = wing_correlations(sig,imgs,fly)
#now do the same thing for the stripe-fixation epoch at the end of the trial
end_cl_idx = flb.idx_by_thresh(stim_cond<0)[-1]
cl_selection_mask = np.zeros_like(flight_mask)
cl_selection_mask[end_cl_idx] = 1
cl_idxs = squeeze(np.argwhere(cl_selection_mask & flight_mask))
imgs = movie[cl_idxs,:,:]
sig = rwing[cl_idxs]
#again correlate with the right wingstroke data
cl_stim_metadata = wing_correlations(sig,imgs,fly)
#now cycle-average the image stream and left and right wing signals
rdict = cycle_average(movie,cycle_idxs,times,flight_mask)
ave_movie = np.mean(rdict['sig_resamp'],axis = 0)
rdict = cycle_average(rwing,cycle_idxs,times,flight_mask)
ave_rwing = np.mean(rdict['sig_resamp'],axis = 0)
rdict = cycle_average(lwing,cycle_idxs,times,flight_mask)
ave_lwing = np.mean(rdict['sig_resamp'],axis = 0)
rdict = cycle_average(ypos,cycle_idxs,times,flight_mask)
ave_ypos = np.mean(rdict['sig_resamp'],axis = 0)
return {'sin_stim_metadata':sin_stim_metadata,
'cl_stim_metadata':cl_stim_metadata,
'ave_movie':ave_movie,
'ave_rwing':ave_rwing,
'ave_lwing':ave_lwing,
'ave_ypos':ave_ypos}
#calculates and save the 'metadata' for each fly.
#import cPickle
#for swarm in [U82_X_J121, U82_X_J18, U82_X_J100,U82_X_J22, U82_X_J23, U56_X_J133]:
# #for swarm in [U82_X_J18]:
# for fly in swarm.flies:
# try:
# #print(fly.fly_num)
# rdict = calc_fly_stats(fly)
# fname = fly.fly_path + 'meta_data.cpkl'
# f = open(fname,'wb')
# cPickle.dump(rdict,f)
# f.close()
# del(rdict)
# except Exception as inst:
# print inst
In [88]:
for swarm_key in swarm_dict.keys():
print swarm_key
group_movie = np.mean(img_dict[swarm_key]['movie_list'],axis = 0)
movie_std = np.std(group_movie[~isnan(group_movie)])
movie_min = 1-movie_std*8
movie_max = 1+movie_std*8
rescaled = (group_movie-movie_min)*255/(movie_std*16)
rescaled[isnan(rescaled)] = 255/2
rescaled[rescaled>255] = 255
rescaled[rescaled<0] = 0
tifffile.imsave(swarm_key + '.tiff',np.uint8(rescaled))
In [86]:
In [83]:
In [186]:
In [186]:
In [186]:
In [98]:
for i,p in enumerate([0,30,60,90]):
subplot(1,4,i+1)
imshow(group_movie[p],cmap = cm.gray,vmin = movie_min,vmax = movie_max)
In [79]:
In [75]:
mov = np.array(img_dict[swarm_key]['movie_list'])
In [103]:
tifffile.imsave('test2.tiff',uint8(mov[1,:,:,:]*50))
In [103]:
imshow(hg3_mask & ~(hg4_mask | i2_mask | hg1_mask))
Out[103]:
In [113]:
shape(img_dict[swarm_key]['rwing_list'])
Out[113]:
In []:
In []:
In []:
In []:
In []:
In []:
In []:
In [19]:
In [20]:
In [21]:
Out[21]:
In [32]:
display(gcf())
In [34]:
plot(rwing[25000])
Out[34]:
In [23]:
fm = np.array(expmnt.exp_record['flight_mask'])
In [37]:
argwhere(~fm)[0]
Out[37]:
In [43]:
imshow(expmnt.exp_record['tiff_data']['images'][int(argwhere(~fm))])
In [46]:
np.mean(np.array(expmnt.exp_record['tiff_data']['images'][int(argwhere(~fm))]),axis = 0)
In [24]:
F = np.mean(np.array(expmnt.exp_record['tiff_data']['images'])[~fm,:,:],axis = 0)
In [25]:
DF = np.array(expmnt.exp_record['tiff_data']['images'])[fm,:,:]
In [26]:
normed = (DF-F)/F
In [34]:
normed[normed<0] = 0
In [37]:
tifffile.imsave('nstream.tiff', uint8(normed[:500]*50))
In [39]:
plot(rwing[fm])
Out[39]:
In [33]:
hist(uint8(normed[:500].ravel()*255))
Out[33]:
In [165]:
swarm = U82_X_J18
fly = swarm.flies[0]
In [169]:
fname = fly.fly_path + 'meta_data.cpkl'
f = open(fname,'rb');data = cPickle.load(f);f.close()
In [170]:
data.keys()
Out[170]:
In [172]:
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname);img_data = cPickle.load(f);f.close()
test_basis = Basis()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
#compose back to Ca image size
A = vstack((A,[0,0,1]))
S = np.array([[0.25,0,0],[0,0.25,0],[0,0,1]])
Ap = dot(S,A)[:2,:3]
In [174]:
imshow(data['ave_movie'][0,:,:])
Out[174]:
In [182]:
In [186]:
In [190]:
imshow(warped_ave_movie[50,:,:])
Out[190]:
In [163]:
In [163]:
In [163]:
In [163]:
In [164]:
imshow(f(np.linspace(0,2,5000)))
In [321]:
tt
Out[321]:
In [123]:
In [123]:
In [123]:
In [196]:
In [197]:
Out[197]:
In [198]:
In [199]:
Out[199]:
In [201]:
Out[201]:
In [123]:
In []:
In []:
In []:
In []:
In []:
To perform ICA I take the first fly in the squadron and warp the movie into a scaled-down version of the confocal frame. The projection of the segemented muscles from the confocal stacks can then be used as pixel masks on the transformed image stack. An alternative would be to warp the masks into the Ca++ imaging frame; the approach I am taking now is conceptually simpler though no doubt computationally much slower. In the analysis these masks are used in two ways: 1) reduce the size of the problem 2) provide spatial information
In []:
First I performed ICA on a mask defined as the union of all the muscles. The 12 columns of the mixing matrx are reshaped into annatomical space and shown as a heat map on the left. The corresponding component is shown to the right. Time is in sec.
In [77]:
ica_filters, ica_sig = st_ica(mask_signal, ncomp = 7, mu = 0.8)
In [89]:
for i in range(shape(ica_filters)[0]):
tmpimg = zeros((256,256))
tmpimg[mask] = ica_filters[i,:]+abs(np.min(ica_filters[i,:]))
figure()
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
thorax_view.plot(little_frame,lw = 0.3)
imshow(tmpimg)
In [16]:
fly =U_X_J.flies[0]
#load the warped movie
warped_movie = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()
#only did one experiment on these flies
expmnt = fly.experiments.values()[0]
#Map some variables
stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])
lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])
rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])
lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])
aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])
xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])
flight_mask = aux < 0.2
trial_selection_mask = stim_cond>0
cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask)
In [21]:
nc = 7
posterior_muscle_mask = i2_mask | iii3_mask \
| iii24_mask | hg1_mask | hg2_mask \
| hg3_mask | hg4_mask
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,posterior_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
posterior_muscle_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
plot_comp_img(A_,posterior_muscle_mask)
In [15]:
def runICA(warped_movie,mask,n_components = 2):
"""get the ICA signals,mixing matrix and unmixing matrix given movie,
mask and expected number of components"""
mask = cv2.pyrDown(cv2.pyrDown(np.float64(mask))) >0
mask_signal = warped_movie[:,mask]
mask_signal = np.float64(mask_signal)
from sklearn.decomposition import FastICA, PCA
ica = FastICA(n_components=n_components,max_iter=1000)
S = mask_signal/mask_signal.std(axis=0)
S_ = ica.fit_transform(mask_signal) # Reconstruct signals
A_ = ica.mixing_ # Get estimated mixing matrix
C_ = ica.components_
return S_,A_,C_
def plot_comp_img(A_,mask):
"""just plot a grid of the component images"""
figure(figsize = (8.5,12))
mask = cv2.pyrDown(cv2.pyrDown(mask.astype(float))) > 0
n_components = shape(A_)[1]
n_rows = int(ceil(nc/3.0))
for i in range(n_components):
col = mod(i,3)
row = i/3
subplot2grid((n_rows,3),(row,col))
tst_img = np.zeros((256,256))
tst_img[mask] = A_[:,i] - np.min(A_)
imshow(tst_img,vmin= 0)#,vmax = 100)
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
thorax_view.plot(little_frame,lw = 0.3)
gca().set_xticks([]);gca().set_yticks([])
gca().set_xbound((0,256));gca().set_ybound((0,256))
def plot_ICA(warped_movie,
mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_):
"""plot a summary of the ICA data from the entire experiment using the components
extracted from a small subset"""
n_components = shape(A_)[1]
n_samples = sum(~flight_mask[:30000])
mask = cv2.pyrDown(cv2.pyrDown(mask.astype(float))) > 0
dat = warped_movie[~flight_mask[:30000],:,:][:,mask]
sigs = np.dot(C_,dat.T)
figure(figsize = (8.5,12))
rwing = rad2deg(rwing/5.0)
#plot the rwing amp cycle average
subplot2grid((n_components+2,5),(0,4))
cd = cycle_average(rwing,cycle_idxs,times)
plot(cd['time_resamp'],cd['sig_mean'])
#plot the rwing example trace
subplot2grid((n_components+2,5),(0,1),colspan = 3)
plot(times[~flight_mask[:30000]],rwing[~flight_mask[:30000]])
ax = gca()
ax.set_ybound((50,90))
raw_sig = np.mean(dat,axis=1)
#plot the cycle average of the raw signal
subplot2grid((n_components+2,5),(1,4))
cd = cycle_average(raw_sig,cycle_idxs,times)
plot(cd['time_resamp'],cd['sig_mean'])
gca().set_yticklabels([])
#plot an example trace of the raw signal
subplot2grid((n_components+2,5),(1,1),colspan = 3,sharex = ax)
plot(times[~flight_mask[:30000]],raw_sig)
gca().set_yticklabels([])
ax = gca()
#itterate through the components
for i in range(n_components):
col = 0
row = i
#first an image of the components
subplot2grid((n_components+2,5),(row+2,col))
tst_img = np.zeros((256,256))
tst_img[mask] = A_[:,i] - np.min(A_)
imshow(tst_img,vmin= 0)#,vmax = 100)
little_frame = copy.copy(thorax.basis)
little_frame['p'] = thorax.basis['p']/4
little_frame['a1'] = thorax.basis['a1']/4
little_frame['a2'] = thorax.basis['a2']/4
thorax_view.plot(little_frame,lw = 0.3)
gca().set_xticks([]); gca().set_yticks([])
gca().set_xbound((0,256)); gca().set_ybound((0,256))
#now the cycle average
subplot2grid((n_components+2,5),(row+2,col+4),colspan = 1)
cd = cycle_average(sigs[i,:],cycle_idxs,times)
plot(cd['time_resamp'],cd['sig_mean'])
gca().set_yticklabels([])
#now an example trace
subplot2grid((n_components+2,5),(row+2,col+1),colspan = 3,sharex = ax)
ax = gca()
plot(times[:n_samples],sigs[i,:])
gca().set_yticklabels([])
ax.set_xbound((40,80))
In [97]:
In [97]:
In [22]:
In [97]:
In [82]:
close('all')
In [97]:
figure();plot(ica_sig[1,:]);
In [29]:
In [30]:
Out[30]:
In []:
In []:
In []:
In []:
In []:
In [26]:
import numpy as np
from numpy import dot,argsort,diag,where,dstack,zeros, sqrt,float,real,linalg
#from numpy import *
from numpy.linalg import eigh, inv, eig, norm
from numpy.random import randn, rand
import time
try:
from scipy.stats import skew
_skew_loaded = True
from scipy.linalg import orth
_orth_loaded = True
except:
_orth_loaded = False
_skew_loaded = False
## TODOs:
## [ ] show how much variance is accounted by pcs
## [X] function to reshape back to movie
## [X] behave nicely if no scipy available (can't sort by skewness then)
## [ ] simple GUI (traits?)
## [ ] masks from ICs or PCs
def pca_svd(X):
ndata, ndim = X.shape
Y = X - X.mean(axis=-1)[:, np.newaxis] # remove mean
U,s,Vh = np.linalg.svd(X, full_matrices=False)
Vh = Vh[:ndata] #only makes sense to return the first num_data
return U,Vh,s
def pca1 (X, verbose=False):
"""
Simple principal component analysis
X as Npix by Nt matrix
X should be normalized and centered beforehand
--
returns:
- EV (Nt by Nesq:esq>0), matrix of PC 'signals'
(eigenvalues of temporal covariance matrix). Signals are in columns
- esq, vector of eigenvalues
"""
tick = time.clock()
Nx, Nt = X.shape
Y = X - X.mean(axis=-1)[:, np.newaxis] # remove mean
C = dot(Y, Y.T)/Nt # temporal covariance matrix
es, EV = eigh(C) # eigenvalues, eigenvectors
## take non-negative eigenvalues
non_neg, = where(es>=0)
neg = where(es<0)
if len(neg)>0:
if verbose:
print "pca1: Warning, C have %d negative eigenvalues" %len(neg)
es = es[non_neg]
EV = EV[:,non_neg]
#S = diag(np.sqrt(es))
#whitenmat = dot(EV, inv(S)) # whitenting matrix
ki = argsort(es)[::-1] # indices for sorted eigenvalues
#U = dot(X, whitenmat) # spatial filters
if verbose:
print "PCA finished in %03.2f sec" %(time.clock() - tick)
return EV[:,ki], es[ki]
# note: returns EVs stored in columns, e.g.
# [[x1, x2, ... xn],
# [y1, y2, ... yn],
# ...]
## Whitening
def whiten_mat1(X):
EV,es = pca1(X)
S = diag(np.sqrt(es))
return dot(EV, inv(S))
def st_ica(X, ncomp = 20, mu = 0.3):
"""Spatiotemporal ICA for sequences of images
Input:
X -- list of 2D arrays or 3D array with first axis = time
mu -- spatial vs temporal, mu = 0 -> spatial; mu = 1 -> temporal
"""
data = X
#data = reshape_from_movie(X) # nframes x npixels
sh = X[0].shape
pc_filters, pc_signals, ev = whiten_svd(data, ncomp)
mux = sptemp_concat(pc_filters, pc_signals, mu)
_, W = fastica(mux, whiten=False)
ica_sig = dot(W, pc_signals)
ica_filters = dot(dot(diag(1.0/np.sqrt(ev)), W), pc_filters)
if _skew_loaded:
skewsorted = argsort(skew(ica_sig, axis=1))[::-1]
else:
skewsorted = range(nIC)
return ica_filters[skewsorted], ica_sig[skewsorted]
def transp(m):
"conjugate transpose"
return m.conjugate().transpose()
pow3nonlin = {'g':lambda X: X**3,
'gprime': lambda X: 3*X**2}
pow3nonlinx = {'g':lambda X,args: X**3,
'gprime': lambda X,args: 3*X**2}
def _sym_decorrelate(X):
"W <- W \cdot (W^T \cdot W)^{-1/2}"
a = dot(X, transp(X))
ev, EV = linalg.eigh(a)
return dot(dot(dot(EV, np.diag(1.0/np.sqrt(ev))),
EV.T),
X)
def _ica_symm(X, nIC=None, guess=None,
nonlinfn = pow3nonlin,
termtol = 5e-7, max_iter = 2e3,
verbose=False):
"Simplistic ICA with FastICA algorithm"
nPC, siglen = map(float, X.shape)
nIC = nIC or nPC
guess = guess or np.random.normal(size=(nPC,nPC))
guess = _sym_decorrelate(guess)
B, Bprev = zeros(guess.shape, np.float64), guess
iters, errx = 0,termtol+1
g, gp = pow3nonlin['g'], pow3nonlin['gprime']
while (iters < max_iter) and (errx > termtol):
bdotx = dot(Bprev, X)
gwtx = g(bdotx)
gp_wtx = gp(bdotx)/siglen
B = dot(gwtx, X.T)/siglen - dot(np.diag(gp_wtx.mean(axis=1)), Bprev)
B = _sym_decorrelate(B)
errx = max(abs(abs(np.diag(dot(B, Bprev.T)))-1))
Bprev = np.copy(B)
iters += 1
if verbose:
if iters < max_iter:
print "Success: ICA Converged in %d tries" %iters
else:
print "Fail: reached maximum number of iterations %d reached"%maxiters
return B.real
def whiten_svd(X, ncomp=None):
n,p = map(float, X.shape)
Y = X - X.mean(axis=-1)[:, np.newaxis]
u, s, vh = linalg.svd(X, full_matrices=False)
K = (u/s).T[:ncomp]
return np.dot(K, X), K, s[:ncomp]
def fastica(X, ncomp=None, whiten = True,
algorithm = 'symmetric',
tol = 1e-3, max_iter = 1e3, guess = None):
n,p = map(float, X.shape)
if whiten:
XW, K, _ = whiten_svd(X, ncomp) # whitened and projected data
else:
XW = X.copy()
XW *= np.sqrt(p)
kwargs = {'termtol':tol, 'nonlinfn':pow3nonlin,
'max_iter':max_iter, 'guess':guess}
algorithms = {'symmetric':_ica_symm, 'deflation':None}
fun = algorithms.get(algorithm, 'symmetric')
W = fun(XW, **kwargs)
if whiten:
S = dot(dot(W,K),X)
else:
S = dot(W,X)
return S, W
def fastica_defl(X, nIC=None, guess=None,
nonlinfn = pow3nonlin,
termtol = 5e-7, maxiters = 2e3):
tick = time.clock()
nPC, siglen = X.shape
nIC = nIC or nPC-1
guess = guess or randn(nPC,nIC)
if _orth_loaded:
guess = orth(guess)
B, Bprev = zeros(guess.shape, np.float64), zeros(guess.shape, np.float64)
iters, minAbsCos = 0,0
errvec = []
icc = 0
while icc < nIC:
w = randn(nPC,1) - 0.5
w -= dot(dot(B, transp(B)), w)
w /= norm(w)
wprev = zeros(w.shape)
for i in xrange(long(maxiters) +1):
w -= dot(dot(B, transp(B)), w)
w /= norm(w)
#wprev = w.copy()
if (norm(w-wprev) < termtol) or (norm(w + wprev) < termtol):
B[:,icc] = transp(w)
icc += 1
break
wprev = w.copy()
return B.real, errvec
### Some general utility functions:
### --------------------------------
def gauss_kern(xsize, ysize=None):
""" Returns a normalized 2D gauss kernel for convolutions """
xsize = int(xsize)
ysize = ysize and int(ysize) or xsize
x, y = mgrid[-xsize:xsize+1, -ysize:ysize+1]
g = np.exp(-(x**2/float(xsize) + y**2/float(ysize)))
return g / g.sum()
def gauss_blur(X,size=1.5):
return signal.convolve2d(X,gauss_kern(size),'same')
def reshape_from_movie(mov):
l,w,h = mov.shape
return mov.reshape(l,w*h)
def reshape_to_movie(X,nrows,ncols):
Np, Nt = X.shape
return X.T.reshape(Nt,nrows,ncols)
def sptemp_concat(filters, signals, mu):
if mu == 0:
out= filters # spatial only
elif mu == 1:
out= signals # temporal only
else:
out = np.concatenate(((1-mu)*filters, mu*signals),
axis = 1)
return out / np.sqrt(1-2*mu+2*mu**2)
def DFoF(X):
"""
Delta F / mean F normalization for each pixel
"""
xmean = X.mean(axis=1).reshape(-1,1)
xmean_z = where(xmean==0)
xmean[xmean_z] = 1.0
out = X/xmean - 1
out[xmean_z,:] = 0
return out
def DFoSD(X, tslice = None):
"""
Delta F / S.D.(F) normalization
"""
if not isinstance(tslice, slice):
tslice = tslice and slice(*tslice) or slice(None)
xsd = X[:,tslice].std(axis=1).reshape(-1,1)
xmean = X[:,tslice].mean(axis=1).reshape(-1,1)
return (X-xmean)/xsd
def demean(X):
"""
Remove mean over time from each pixel
Frames are flattened and are in columns
"""
xmean = X.mean(axis=0).reshape(1,-1)
return X - xmean
def winvhalf(X):
"raise matrix to power -1/2"
e, V = eigh(X)
return dot(V, dot(diag((e+0j)**-0.5),transp(V)))
def msqrt(X):
e, V = eigh(X)
return dot(V, dot(diag(((e+0j)**0.5)), transp(V)))
def mpower(M, p):
"""
Matrix exponentiation, works for Hermitian matrices
"""
e,EV = linalg.eigh(M)
return dot(transp(EV),
dot(diag((e+0j)**p), EV))
def mask4overlay(mask,colorind=0):
"""
Put a binary mask in some color channel
and make regions where the mask is False transparent
"""
sh = mask.shape
z = np.zeros(sh)
stack = dstack((z,z,z,np.ones(sh)*mask))
stack[:,:,colorind] = mask
return stack
def flcompose2(f1,f2):
"Compose two functions from left to right"
def _(*args,**kwargs):
return f2(f1(*args,**kwargs))
return _
def flcompose(*funcs):
"Compose a list of functions from left to right"
return reduce(flcompose2, funcs)
In []:
In []:
In []:
In []:
In []:
In []:
In []:
In [184]:
nc = 3
#ica_data = warped_movie[~flight_mask[:10000],:,:]
S_,A_,C_ = runICA(ica_data,b3_mask,n_components = nc)
plot_ICA(warped_movie,
b3_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
plot_comp_img(A_,b3_mask)
In [190]:
print mean(A_,axis = 0)
In [165]:
In [165]:
In [165]:
In [165]:
In []:
In [19]:
nc = 12
all_muscle_mask = b1_mask | b2_mask \
| b3_mask | i1_mask \
| i2_mask | iii1_mask \
| iii3_mask | iii24_mask \
| hg1_mask | hg2_mask \
| hg3_mask | hg4_mask
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,all_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
all_muscle_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
plot_comp_img(A_,all_muscle_mask)
In [31]:
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
| b3_mask | i1_mask \
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,anterior_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
anterior_muscle_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
figure()
plot_comp_img(A_,anterior_muscle_mask)
In [32]:
nc = 2
ica_data = warped_movie[~flight_mask[:5000],:,:]
b3_bool_mask = b3_mask & ~(i1_mask | b1_mask | b2_mask)
S_,A_,C_ = runICA(warped_movie,b3_bool_mask,n_components = 2)
plot_ICA(warped_movie,
b3_bool_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
#display(gcf());close()
In [44]:
fly = U_X_J.flies[0]
#load the warped movie
warped_movie = tifffile.TiffFile(fly.fly_path + 'warped_imgs.tiff').asarray()
#only did one experiment on these flies
expmnt = fly.experiments.values()[0]
#Map some variables
stim_cond = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['DAC2'])
lwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph0'])
rwing = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph1'])
lmr = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph2'])
aux = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ph3'])
xpos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Xpos'])
ypos = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['Ypos'])
times = np.array(expmnt.exp_record['tiff_data']['axon_framebase']['times'])
flight_mask = aux<0.0
trial_selection_mask = stim_cond>0
cycle_idxs = collect_cycles(ypos,stim_cond,trial_selection_mask)
In [45]:
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
| b3_mask | i1_mask \
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,anterior_muscle_mask,n_components = nc)
plot_ICA(warped_movie,
anterior_muscle_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
figure()
#plot_comp_img(A_,anterior_muscle_mask)
Out[45]:
In [46]:
nc = 3
#i1_mask = b1_mask | b2_mask \
# | b3_mask | i1_mask \
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,i1_mask,n_components = nc)
plot_ICA(warped_movie,
i1_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
figure()
#plot_comp_img(A_,anterior_muscle_mask)
Out[46]:
In [48]:
nc = 3
#i1_mask = b1_mask | b2_mask \
# | b3_mask | i1_mask \
ica_data = warped_movie[~flight_mask[:5000],:,:]
S_,A_,C_ = runICA(ica_data,b2_mask,n_components = nc)
plot_ICA(warped_movie,
b2_mask,
flight_mask,
times,
rwing,
cycle_idxs,
S_,A_,C_)
figure()
#plot_comp_img(A_,anterior_muscle_mask)
Out[48]:
In [47]:
draw()
In [34]:
plot(proj[1,:])
Out[34]:
In [57]:
plot(cycle_average(S_[:,0],cycle_idxs,times)['sig_median'])
Out[57]:
In [126]:
s1 = sum(warped_movie[:,mask],axis = 1)
In [122]:
i1_bool_mask = i1_mask & ~(b2_mask)
mask = cv2.pyrDown(cv2.pyrDown(np.float64(i1_bool_mask))) >0
s2 = sum(warped_movie[:,mask],axis = 1)
In [123]:
from sklearn.decomposition import FastICA
ica = FastICA(n_components=1)
S = c_[s1,s2]/c_[s1,s2].std(axis = 0)
S_ = ica.fit_transform(c_[s2,s1]) # Reconstruct signals
A_ = ica.mixing_ # Get estimated mixing matrix
C_ = ica.components_
In [127]:
imshow(b3_bool_mask)
Out[127]:
In [119]:
plot(flight_mask)
Out[119]:
This approach seems like it is useful in the sense that it identifies neworks of activity, but it is difficult to interpret in terms of isolating the signal from individual muscles. This is likely due to the fact that optimization of the mixing matrix is not constrained by any spatial information. I have found several papers that describe ICA with spatial constraints, but implementing these require that I roll my own ICA optimization algorithm - and this would take some time. As a first pass I took an ad-hoc approach to applying spatial constraints: first by breaking the analysis up into the anterior and posterior cluster of muscles, then looking at the components of each muscle one-by-one.
In [15]:
#the anterior cluster
nc = 4
anterior_muscle_mask = b1_mask | b2_mask \
|b3_mask | i1_mask
S_,A_ = runICA(anterior_muscle_mask,n_components = nc)
plot_ICA(anterior_muscle_mask,S_,A_,n_components = nc);display(gcf());close()
In [17]:
nc = 3
b2_bool_mask = b2_mask
S_,A_ = runICA(b2_bool_mask,n_components = nc)
plot_ICA(b2_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [18]:
nc = 3
i1_bool_mask = i1_mask
S_,A_ = runICA(i1_bool_mask,n_components = nc)
plot_ICA(i1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [19]:
nc = 3
b3_bool_mask = b3_mask
S_,A_ = runICA(b3_bool_mask,n_components = nc)
plot_ICA(b3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [20]:
nc = 3
b1_bool_mask = b1_mask
S_,A_ = runICA(b1_bool_mask,n_components = nc)
plot_ICA(b1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [21]:
nc = 3
iii3_bool_mask = iii3_mask
S_,A_ = runICA(iii3_bool_mask,n_components = nc)
plot_ICA(iii3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [22]:
nc = 3
iii24_bool_mask = iii24_mask
S_,A_ = runICA(iii24_bool_mask,n_components = nc)
plot_ICA(iii24_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [23]:
nc = 3
hg1_bool_mask = hg1_mask
S_,A_ = runICA(hg1_bool_mask,n_components = nc)
plot_ICA(hg1_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [24]:
nc = 3
hg2_bool_mask = hg2_mask
S_,A_ = runICA(hg2_bool_mask,n_components = nc)
plot_ICA(hg2_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [25]:
nc = 3
hg3_bool_mask = hg3_mask
S_,A_ = runICA(hg3_bool_mask,n_components = nc)
plot_ICA(hg3_bool_mask,S_,A_,n_components = nc);display(gcf());close()
In [26]:
nc = 3
hg4_bool_mask = hg4_mask
S_,A_ = runICA(hg4_bool_mask,n_components = nc)
plot_ICA(hg4_bool_mask,S_,A_,n_components = nc);display(gcf());close()
it is also possible to define masks that atempt to excude pixels from overlaping regions
In [27]:
nc = 3
b2_only_boolmask = i1_mask & ~(b2_mask | b1_mask |b3_mask)
S_,A_ = runICA(b2_only_boolmask,n_components = nc)
plot_ICA(b2_only_boolmask,S_,A_,n_components = nc);display(gcf());close()
In [0]:
In [0]:
In [64]:
In [45]:
### function to itterate through a large list of flies
### and make figures of the fits - probably should migrate this to
### an independent script
import db_access as dba
import flylib as flb
import cv2
import cPickle
fly_db = dba.get_db()
def plot_maxfits(flylist):
numfigs = np.ceil(len(flylist)/6.)
chunks = np.array_split(flylist,numfigs)
for chunk in chunks:
fit_swarm = flb.Squadron(fly_db,chunk)
figure(figsize = (8.5,11))
for i,fly in enumerate(fit_swarm.flies):
output_shape = (1024,1024)
#fly = swarm.flies[0]
tiffdata = np.array(fly.fly_record['experiments'].values()[0]['tiff_data']['images'])
pkname = fly.fly_path + '/basis_fits.cpkl'
f = open(pkname)
img_data = cPickle.load(f)
f.close()
test_basis = Basis()
for key in ['A','p','a1','a2']:
test_basis[key] = img_data[key]
src_p0 = test_basis['a1'] + test_basis['p']
src_p1 = test_basis['p']
src_p2 = test_basis['a2'] + test_basis['p']
dst_p0 = confocal_basis['a1'] + confocal_basis['p']
dst_p1 = confocal_basis['p']
dst_p2 = confocal_basis['a2'] + confocal_basis['p']
maximg = np.max(tiffdata,axis=0)
col = mod(i,2)
row = i/2
subplot2grid((3,2),(row,col))
A = cv2.getAffineTransform(np.float32([src_p0,src_p1,src_p2]),np.float32([dst_p0,dst_p1,dst_p2]))
imshow(cv2.warpAffine(maximg.T,A,output_shape),cmap = cm.gray)
thorax_view.plot(thorax.basis,alpha = 0.5)
gca().set_ybound(0,1024)
gca().set_xbound(0,1024)
gca().set_xticks([])
gca().set_yticks([])
gca().set_title(('Fly%04d')%(fly.fly_num))